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Abstract 


There  are  several  very  fast  direct  methods  which  can  be  used  to 
sol/e  the  discrete  Poisson  equation  on  rectangular  domains.  We  show  that 
these  methods  can  also  be  used  to  treat  problems  on  irregular  regions. 


1.  Introduction.  Within  the  past  few  years,  several  very  fast  and  accurate 


direct  methods  have  been  developed  for  solving  finite  difference  approximations 
to  the  Poisson  equation, 

Au  =  f  in  R  , 
u  =  g  on  Sr  .  J 

These  methods  can  usually  be  applied  only  on  rectangular  regions,  although  the 
differential  operator  and  boundary  conditions  can  be  more  general  than  those  in 
the  Poisson  equation.  In  this  paper,  we  will  show  how  these  algorithms  for 
rectangular  domains  can  also  be  used  effectively  on  irregular  regions .  The 
approach  used  is  similar  to  that  employed  by  Hockney  [l6,  17],  Buneman  [7], 
and  George  [lU].  We  also  mention  the  work  of  Angel  [1-^],  Angel  and  Kalaba  [^], 
Collins  and  Angel  [9],  Kalaba  [20],  end  Roache  [22]  on  the  use  of  direct  methods 
for  problems  in  irregular  regions. 

We  will  not  discuss  the  details  of  any  specific  direct  method.  A  survey 
of  these  procedures  is  given  in  [11],  and  in  particular  we  cite  the  recent  work 
of  Buneman  [6],  Buzbee,  Golub,  and  Nielson  [8],  and  Hockney  [16]. 

We  will  also  not  consider  the  derivation  of  the  finite  difference  equations 
that  approximate  the  partial  differential  equation.  This  subject  is  treated  in 
detail  by  Forsythe  6ind  Wasow  [15]^  and  we  assume  that  the  problem  has  been 
reduced  to  finding  the  solution  of  a  matrix  equation  Ax  =  y  .  The  matrix  A 
is  frequently  very  large  and  sparse,  but  its  structure  does  not  permit  the 
application  of  the  most  efficient  direct  methods.  For  our  computational  procedure, 
we  alter  certain  rows  of  A  to  obtain  a  matrix  B  ,  and  we  will  show  how  to 
define  a  modified  right-hand  side  z  so  that  the  solution  x  also  satisfies  the 
equation  Ex  =  z  .  The  matrix  B  is  chosen  so  that  these  equations  can  be 
solved  by  the  direct  methods. 
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This  method  is  computationally  advantageous  when  we  are  solving  a 
sev^uence  of  equations  Ax^  =  y.  .  This  situation  frequently  arises  in 
time-dependent  partial  differential  equations,  in  nonlinear  problems,  and  in 
linear  problems  where  the  right-hand  side  is  varied  but  the  region  and 
differential  operator  remain  the  san.e.  After  some  initial  computation,  each 
solution  X.  can  be  obtained  in  approximately  twice  tne  time  required  for 
the  solution  of  an  equation  Bx  =  z 

In  Gectlons  2  and  3  we  derive  this  algorithm  in  a  general  foirni.  We 
describe  a  number  of  applications  of  the  method  in  Sections  k  and  5,  and  in 
Section  t  we  present  some  computational  results. 
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2.  Method  of  Solution  if  det  B  ^  0  .  Suppose  that  we  are  given  an  n  hy  n 
matrix  A  and  an  integer  p  with  1  <  p  <  n .  We  wish  to  modify  p  rows  of  A 
to  obtain  another  matrix  B  .  Without  loss  of  generality  we  sissume  that  the 
first  p  rows  of  A  are  to  he  changed,  since  we  can  achieve  this  situation 
hy  multiplying  A  hy  a  suitable  permutation  matrix.  However,  we  enqphasize 
that  this  multiplication  Should  not  he  done  explicitly  in  the  coeiputational 
procedure.  Rather,  the  rearrangement  of  rows  Should  he  done  implicitly  hy 
indexing.  IHie  direct  methods  mentioned  later  in  the  paper  require  that  B 
has  a  particular  structure,  which  could  he  altered  hy  the  permutation 
transformation. 

Partition  A  in  the  form 


where  A^  is  a  p  hy  n  matrix  and  is  an  (n-p)  hy  n  matrix.  We  then 
write 


where  B^  is  a  p  hy  n  matrix.  For  the  remainder  of  this  section,  we 
sussume  that  det  B  /  0  . 

Sxqipose  we  are  given  a  linear  equation  Ax  =  y  .  We  partition  y  in 
the  same  way  as  A  ,  and  write 
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Let  y  be  any  vector  of  the  form 


If  W  is  an  arbitrary  nonsingular  p  by  p  matrix,  we  define  an  n  by  p 
matrix  W  by 


Define  the  p  by  p  matrix  C  by 


C  =  M  . 


Following  Hockney  [l6],  we  call  C  the  capacitance  matrlXT^  Assume  that  there 
exists  a  p  by  1  vector  ^  that  is  a  solution  to  the  equation 


CP-y  . 

^  X  ^ 


(1) 


Since  A  6md  B  differ  only  in  the  first  p  rows,  it  is  easy  to  verify  that 
a  solution  x  to  the  equation  Ax  =  y  is  given  by 


X  =  (y  +  WP)  . 


We  first  show  that  this  method  of  obtaining  the  solution  x  will  be  valid 
whenever  the  original  system  Ax  =  y  is  consistent. 

Theorem.  ^  det  B  /  0  ,  then 


detC  =■  ('irtA)(drt1l)  _ 

D 


*1  -1 

-  Hockney  actually  refers  to  C  as  the  capacitance  matrix.  Since  C  may 
be  singular  in  o\ir  development,  we  have  adopted  the  present  notation. 


If  thft  sy>te&  Ax  B  y  is  consistent,  then  Eq.  (l)  Is  also  consistent. 

^  ^  ■—  ■  ■  "  --  BM  >■•■>  1  ' 

Proof.  Partition  in  the  fora 


B 


"■(»X  ■>.). 


ahere  is  n  by  p  and  is  n  by  (n-p)  .  It  then  follows  that 


BB 


-1 


/I  0 


and 


AB 


-1 


.®1 


‘^2®1  ^2®2^  ^  ° 


A^Dl  A^Lg, 
I 


Thus  we  have 


C  =  Aj^B'^W  «  • 


and  so 


det  C  -  det  (A^  D^)  det  V 

«  det  (AB*^)  detW 

(det  A)  (det  W) 

‘  detB 


To  prove  the  consistency  statenenty  sv^ipoee 


Write 
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axid  define  an  n  by  1  vector  x  ^7 


Since  the  system  Ax  s  y  is  assuBied  to  be  consistent,  we  therefore  have 
*T  T 

r  y  =  0  ,  which  is  the  sane  as  y  -  ^^2^  -  0  .  But  then 

l{h.  ■  '  I*(2i  -  -  ®iZs) 

-  0  , 

Which  is  the  consistency  condition  for  Eq.  (1). 

The  Woodbury  formila  [18,  pp.  123  -124]  for  the  inverse  of  a  natrlx 
(B  +  F  G)  is 

^B  +  Fg)  =  B"^  -  F  (l  +  GB“^f)  ^  GB“^) 

This  equation  has  been  used  in  direct  nethods  for  solving  the  Poisson  equation 
by  George  [l4],  and  for  the  blhannonlc  equation  by  Golifl)  [13].  If  A  is  non¬ 
singular  we  write 

A  =  B  +  FG  , 
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ithere  ¥  *  ¥  ,  and  6  is  the  p  hy  n  natrix  given  by 

G  -  (A^  -  Bj)  . 

For  the  case  in  lAdch  A  is  nonsingalar,  the  algorithm  we  have  derived  is 
equivalent  to  using  the  Woodbury  fonmila  for  A"^ . 

Suppose  that  we  have  a  very  efficient  method  for  solving  equations  of 
the  form  Bz  »  w  .  The  solution  of  the  equation  Ax  s  y  then  proceeds  in 
the  following  steps: 

(1)  Conpute  C  -  A^B'^W  , 

(2)  Cospute  X  »  B"^y  , 

(3)  Solve  the  equation  C  S  »  y.  -  A.  x 

The  solution  x  can  then  be  obtained  from  the  fonmila 

see 

x-B’^(y  +  Wp)  .  (2) 

If  it  is  possible  to  store  the  vector  x  and  the  matrix  B  s  B'^W  .  then  x 
can  also  be  coBQnzbed  trom 

X  .  X  +  BP  .  (3) 

The  decision  Whether  to  use  Eq.  (2)  or  Eq.  (3)  would  be  made  on  consideration 
of  storage  requirements,  and  on  the  relative  sp-ted  of  solving  the  system  in 
Eq.  (2)  versus  multiplying  by  the  matrix  in  Eq.  (3).  For  problems  arising 
treat  elliptic  difference  equations,  it  is  frequently  better  to  use  Eq.  (2) 
because  B  has  a  band  structure,  but  the  matrix  B  may  be  full. 

The  type  of  application  we  have  in  mind  for  this  method  is  one  in  Which 
we  have  to  solve  a  mmber  of  equations  A^  ~  ^  case,  we  conpute 

the  capacitance  matrix  and  factor  it  m  part  of  a  preprocessing  stage.  The 
solution  of  eaOh  equation  A^  ^  then  approximately  as  fast  as  the  time 
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it  takes  to  solve  two  equations  Bz  =  v  . 

To  be  specific,  let  6(n)  denote  the  number  of  arithmetic  operations 
necessary  to  solve  a  system  B  z  =  w  .  Then  to  cocipate  C  and  form  its  LU 
decomposition  in  a  preprocessing  stage  requires  approximately 

2  5 

pe(n)  +  k^p  n  +  kgp-' 

operations  (cf.  [19,  Sec.  2.l]).  In  many  cases  the  matrix  is  sparse, 
and  this  estimate  is 

pe(n)  +  kjp^  (4) 

operations.  To  coeqnite  the  solution  to  a  particular  equation  Ax  >  y  using 
£q.  (2)  takes  an  additional 

o 

2e(n)  +  kj^pn  +  k^p 

operations.  If  A^  is  sparse  and  ve  let  W  >  I  ,  this  estimate  can  be  replaced 
by 

2  0(n)  +  kgp^  (5) 

operations.  To  coopute  a  particular  solution  using  £q.  (3)  requires 

p 

0(n)  +  k^pn  +  kgp  (6) 

operations.  In  general  this  estimate  cannot  be  reduced,  because  the  matrix 
B  may  be  full. 
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3.  Method  of  Solution  If  rank(B)  =  n-1  »  The  method  derived  in  Section  2 
gives  a  procedure  for  finding  a  p  by  1  vector  &  such  that  a  solution  x 
to  Ax  =  y  also  satisfies  the  equation 

Bx  =  y  + 

If  B  is  singular,  it  may  not  be  possible  to  find  such  a  vector  &  .  To  show 
T 

this,  suppose  B  v  =  0  but  v  /  0  .  In  order  for  &  to  exist,  we  must 
satisfy  the  consistency  condition 

liy  +  fi^  "  ° 

T  T 

If  V  e^  =  0  for  1  <  i  <  p  and  v  y  /  0  ,  it  is  not  possible  to  satisfy 

Eq.  (7)  •  However,  if  A  is  nonsingular  this  difficulty  does  not  arise, 

T  T 

because  then  the  only  vector  v  satisfying  B  v  =  0  and  v  e^  =  0  for 
l<i<p  is  v  =  0. 

We  will  now  describe  an  algorithm  we  have  used  when  rank(B)  =  n-1  and 
A  is  nonsingular.  There  are  two  advantages  in  treating  this  particular  case. 
First,  the  construction  is  quite  simple,  and  it  is  easy  to  see  how  the  method 
could  be  extended  to  a  more  general  matrix  B  .  Second,  the  case  rank(B)  =  n-1 
has  a  special  significance  in  the  solution  of  partial  differential  equations, 
because  this  condition  is  satisfied  by  the  matrix  corresponding  to  the 
Neumann  problem.  For  simplicity,  we  assume  that  the  matrix  W  of  Section  2 
is  the  identity  matrix. 

Theorem  2.  Assume  that  A  i£  nonsingular  and  rank(B)  =  n-1  ,  aM  let  u 

T 

and  V  ^  two  non-zero  vectors  satisfying  Bu  =  B  v  =  0  .  Then  there  exists  an 

T 

integer  k  with  1  <  k  <  p  such  that  v  e  /  0  .  Define  a  constant 
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and  let  x  be  a  solution  to 

solution  to 

} 

p  matrix  whose  i-th  column  is 
the  vector  •  Then  C  ^  nonsingular ^  and,  if  3  is  the  solution  to 

Cp  =  Yi  -  ^1  X  , 

the  solution  x  Ax  =  y  ^  given  ^ 

P 

X  =  X  +  r  3  Ti  . 

~  "  ^^l  ^ 

Proof.  If  we  partition  v  in  the  same  way  as  y  ,  we  have 

rn  m  fp  T  T  T  T* 

A  V  =  A^  v^  +  A^  v^  and  B  v  =  v^^  +  A^  v^  .  Tnus  if  B  v  =  0  and 

T 

v^  =  0  we  would  have  A  v  =  0  .  Since  A  is  nonsingular  and  y  f  0  this 
cannot  happen,  and  hence  y^  f  0  . 

To  prove  that  C  is  nonsingular,  we  show  that  Cp  =  0  implies  p  =  0  . 

P 

Suppose  p  is  an  arbitrary  vector  such  that  Cp  =  0  .  Then  x  =  Pi  Di 
satisfies  Ax  =  0  ,  and  hence  x  =  0  .  This  implies  that  Bx  =  0  ,  or 

-i)-"  ■ 


Bx  =  y  -  (a  V  y)  e^^  • 

For  1  <  i  <  p  and  i  /  k  let  be  a 

m 

BTl.  =  e.  -  (a  v'^  e.)  e, 
and  let  Tl  =  u  .  Let  C  ^  the  p  ^ 
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Thus  3^  =  0  for  1  <  i  <  P  and  i  /  k  ,  and  the  condition  x  =  0  then 
implies  that  =  0  .  Thus  p  =  0  ,  and  so  C  is  nonsingular. 


Remark.  As  we  discussed  in  Section  2,  the  computation  proceeds  in  the  following 
steps: 

(1)  Compute  (and  factor)  C  , 

(2)  Compute  x  , 

(3)  Solve  for  p  . 

The  solution  x  can  then  be  obtained  from  the  foita’ula 

f-,  h  l  • 


X  =  X  + 


i=l 


However,  if  tie  problem  arises  from  a  partial  differential  equation,  it  is  more 
efficient  computationally  to  obtain  x  in  the  form 


A  A 

X  =  X  +  P  U  , 


where  x  is  a  solution  to 


Bx  =  y  -  (a  v-^  y)  e  +  V  p  (e 

. ^  i=l  ^ 

i/k 


and 


P  =  [u^  J  +  f  p  u^  Tl.  -  u^  x][u^  u]'^  . 
~  ~  1  - 
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4.  Applications  to  Partial 


z.  Suppose  we 


are  given  a  two-dlnensional  bounded  regi.m  R  in  the  x  -  y  plane,  and  we  wish 
to  find  a  solution  u  to  the  Poisson  equation, 

Au  =  f  in  B  A 

u  =  g  on  5r  .  j 

We  assume  that  this  differential  equation  is  approximated  by  a  finite  differ¬ 
ence  equation  (cf.  Forsythe  and  Waaow  [13]).  Thus  we  have  a  finite  set  of 
unknowns  (  |  1  <  i  ^  }  which  approximate  the  solution  u  at  the  grid 

points.  If  we  denote  by  a  finite  difference  approximation  to  the  Laplaclan 
operator  A  ,  by  Rj^  the  discrete  interior  of  the  grid,  and  by  ^  R^^  the  discrete 
boundary  of  the  grid,  then  the  discrete  Poisson  equation  ca:  be  written  in  the 


f  in 


U  =  g  on  SRjj 


:  } 


Let  R^  be  a  discrete  rectangular  region  such  that  c  jRj^  and 
SRj^cR^USr^  ,  and  let  ^  ^  ^  .  Extend  the  functions  f  and  g 

to  the  regions  R^^  and  ^R^  U  dR^  respectively,  and  consider  the  equation 


V' 


.  I 

uK  •  J 


U=g  onSj^UoRj^  .  J 

We  will  solve  Eq.  (9),  and  the  solution  U  will  then  also  satisfy  Eq.  (8). 

Equation  (9)  Is  a  linear  equation  in  the  unknowns  {U^|l^i^n}. 
Observe  that  we  may  have  increased  the  number  of  unknowns  by  the  laibedding 
process,  so  that  n^  x  n  .  We  write  Eq.  (9)  as  a  matrix  equation  AU  =  V  , 


I 


and  the  matrix  A  can  frequently  be  chosen  to  be  block  tr.'.dia^onal  with 
tridiagonal  matrices  as  the  non- zero  blocks  (cf.  [15]). 

Let  p  be  the  number  of  grid  points  in  3^  .  We  modify  the  p  rows  of 
A  and  V  corresponding  to  the  equations 

U  =  g  > 


and  replace  them  with  the  equat  ions 


=  f  on 


This  defines  a  new  matrix  B  and  a  new  ri^t-hand  side  V  .  An  equation 
BU  =  V  corresponds  to  the  difference  equation 


A^U=  f 
U  =  g 


in 

on  S 


(10) 


Since  R^^  is  a  rectangular  region,  we  have  very  fast  methods  for  solving 
Eq.  ( 10) .  We  can  now  apply  the  method  of  Section  2  to  solve  the  equation 
AU  =  V  by  using  the  modified  matrix  B  . 

To  illustrate  this  construction,  let  R  be  a  rectangular  region  with 
an  interior  rectangle  removed,  such  as  that  shown  in  Figure  1,  For  simplic¬ 
ity,  we  assume  that  the  discrete  ioundary  dRj^  is  a  subset  of  SR,  The 
inbedding  rectangle  is  ^  \  U  U  T^^  .  The  only  function  extension  re¬ 
quired  for  this  example  is  that  f  be  defined  (arbitrarily)  in  ^  ^  . 

To  define  this  extension,  we  can  set  f  s  0  in  Sj^  U  T^ ,  or  we  can  define  f 
so  that  it  is  continuous  in  all  of  R^  .  The  advantage  of  using  e.  continuous 
f  is  that  the  solution  to  Eq.  (1  )  is  then  smooth.  However,  the  direct 
methods  used  to  solve  Eq.  (li)  are  so  accurate  that  the  smoothness  of  the 


solution  does  not  appear  to  influence  the  caeqnitational  results.  Therefore, 


In  tho  «xaiiQ>les  we  hare  considered,  we  extend  f  by  setting  f  *  0  in 

If  we  let  W  >  I  in  the  method  of  Section  2,  this  algorithm  is  closely 
connected  with  the  discrete  Green's  function  for  the  region  (cf.  [13  , 
pp.  31^-33JB]).  In  fact,  the  method  is  then  equivalent  to  adding  suitable 
multiples  of  the  discrete  Green's  function  for  the  points  on  so  that 
the  boundary  conditions  on  will  be  satisfied.  Since  we  have  Dirichlet 
boundary  conditions  on  ,  by  a  proper  ordering  of  the  unknowns  we  can  write 

\=  il  0)  . 

Since  B  is  positive  definite  and 

C  =  (I  0)  B*^  /  I 

\  0 

we  see  that  C  is  also  positive  definite  in  this  case.  This  is  advantageous 

T 

l)ecause  CholesI^  decoBopositlon  can  then  be  used  to  ccnpute  an  ll  decomposi- 
tion  of  C  (cf.  [12,  Chap.  23]). 

If  the  grid  on  has  N  points  on  a  side,  we  have  n  =  .  In  that 

case,  we  can  solve  the  system  BU  ■  V  in  approximately 

e(N)  =  5N^  loggN 

operations  (cf.  [U,  p.  260]).  The  preprocessing  then  takes 

5pN^log2N  +  kjp5 

operations  (cf.  Eq.  (4)).  To  solve  Bq.  (8)  for  a  particular  choice  of  f 
and  g  by  using  Eq.  (2)  with  W  »  I  takes  an  additional 

U 


lON^loggN  + 

operations  (cf.  Eq.  (5)).  If  we  use  Eq.  (3)  to  cooipute  the  solution,  it 
takes  an  additional 

5N^log2N  +  k^p/  +  kgp^ 

operations  (cf.  Eq.  (6)).  Thus  if  p  »  loggN  it  is  faster  to  use  Eq.  (2) 
to  compute  the  solution.  We  also  observe  that  for  this  problem  the  matrix 
b“^  is  full  [2:5,  p.  85],  so  to  store  B  in  using  Eq.  (3)  wtjuld  require  pU^ 
locations.  Thus  for  large  values  of  p  and  N  it  is  both  faster  and  more 
economical  in  terms  of  storage  to  use  Eq.  (2)  to  compute  the  solution  to  a 
particular  equation. 

It  shofuld  be  clear  that  the  imbedding  procedure  can  be  applied  to  other 
elliptic  difference  operators  with  other  types  of  boundary  conditions.  To 
be  a  practical  procedure,  we  simply  require  that  we  have  a  fast  method  for 
solving  the  imbedded  problem  in  the  rectmagular  region. 

As  ariother  example,  consider  the  region  shown  in  Figure  2.  This  problem 
arises  in  tne  time-dependent  study  of  a  rotating  fluid  [10],  and  the  fluid  surface 
is  moviiig  slowly.  We  are  given  Dirichlet  boundary  data  on  ,  and  Neumann 
boundary  data  on  .  Tlie  imbedding  rectangle  is  \  U  U  T^^  ,  and 

we  use  Neumann  bour.dar:/  conditions  on  .  Thus  B  corresponds  to  the  Neumann 

problem  on  K^'  ,  and  the  rank  of  B  is  n  -  1  .  The  method  of  Section  3  can  then 
be  applied,  and  direct  methods  for  solving  'che  rectangular  Neumann  problem  are 
given  in  [3]. 

For  an  example  with  the  Poisson  operator  in  another  geometry,  consider 
the  region  in  the  z  -  r  plane  shown  in  Figure  5.  This  problem  arises  in  the 
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time-dependent  study  of  a  plasma  [21i,  and  a  Poisson  equation  must  be  solved 
at  each  time  step.  The  boundary  conditions  are  Dirichlet  on  and  Be-umarn 

on  jjjCl)  .  vie  use  Beumann  boundary  conditions  on  ^  and  “'<1 

Dirichlet  boundary  conditions  on  ^  and  imbedding 

region  B,;  ■=  IS,  U  U  T^  •  The  elliptic  difference  equation  in  R^;  is  solved 

by  the  method  of  matrix  deconposition  [8]. 


5.  Applications  to  Partial  Differential  Equations  "by  Splitting.  There  sire 
many  prohlems  for  which  the  imbedding  approach  is  not  an  economical  algorithm. 
For  example,  imbedding  the  region  in  a  rectangle  may  introduce  an  excessively 

large  number  of  additional  unmiowns  that  are  not  necessary  to  the  solution  of 

the  original  problem.  Another  instance  is  one  in  which  the  differential  oper¬ 
ator  or  the  mesh  size  changes  in  different  parts  of  the  region.  In  this 
section,  we  give  two  such  examples.  In  each  case,  the  method  of  Section  2 
can  he  used  to  split  the  problem  into  two  rectangular  problems,  which  can  be 
solved  by  the  usual  direct  methods. 

Consider  the  elongated  L-shaped  region  in  Figure  4,  and  the  equation 

=  f  in  , 

U  =  g  on  . 

We  assume  that  points  on  the  line  marked  are  all  grid  points.  To  define 

the  matrix  B ,  we  replace  the  equations 

=  f  on 


by  the  equations 


U  =  g  on  Tj^  , 


where  g  has  been  (arbitrarily)  extended  to  Tj^ .  The  solution  of  an  equation 
BU  =  V  now  consists  of  solving  the  two  rectangular  problems 


U.  =  f 

1 

in  b(‘)  , 

U,  =  g 

on  , 

for  i  =  1  ,  2  .  We  can  then  apply  the  method  of  Section  2  to  solve  the  original 


problem.  This  algorithm  is  similar  to  one  developed  in  [  8,  Sec.  9]  for 
non-rect angular  regions. 


As  another  example,  consider  the  multiple-material  problem  shown  in 

Figure  5.  The  differential  equation  is 

S  /  >  V  Su^ 
^  ) 

1  ^  ^  (t(I') 

II 

f  (x  ,  y)  , 

and 

1 

in  R^^^ 

> 

o(x)  » 

1 

in  R^^^ 

• 

The  functions  a^(x),  t  T(y)  are  assumed  to  be  smooth.  Dirichlet 

data  is  given  on  S  R  ,  and  we  require  that  0  be  continuous  euaross  the 
boundary  between  and  R^^^  .  The  coeqnxtational  procedure  is  essentially 

the  same  as  that  for  the  L-shaped  region.  The  only  difference  is  that  in 
forming  the  matrix  B  we  replace  the  equations  for  the  continuity  of 
across  the  line  Tj^  by  the  equations 

U  =  g  0“  • 

As  before,  the  equation  £U  =  V  corresponds  to  the  two  rectangular  problems 

\  (x,  7)  =  b(x  ,  y)  on  Br^^^  , 

for  i  =  1 ,  2  .  These  problems  can  be  solved  directly  by  the  method  of  matrix 
decoiiq)Osition  [  3  ^  See.  8].  A  similar  method  can  be  used  for  the  cause  in 
which  T(y)  is  only  piecewise  smooth. 
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It  1.  clear  that  thla  aplittlng  method  can  be  aEplied  to  the  Poiason 
aquation  ^  regloa.  mich  aa  tlu«  in  Pignre  5  »hen  different  meah  aiaea  are 
uaed  in  8  and  8*  '  .  Ihe  method  derelopcd  in  [  8 ,  Sec.  8]  can  alao  be 
adapted  to  include  rectangular  problema  with  irregular  meahea. 
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6.  Congmtatlonal  Results.  In  Table  1  we  have  tabulated  some  con^nitational 

results  for  two  regions  of  the  form  of  Figure  1.  In  each  case,  a  square 

with  sides  of  length  1  has  a  sjnaaetrically  located  square  removed  from  its 

center.  For  region  1  the  inner  square  has  sides  of  length  and  for 

region  2  the  inner  sides  are  of  length  ,  We  solve  the  Poisson  equation 

2  2 

with  Dirichlet  boundary  conditions  for  the  function  u(x  ,  y)  =  x  +  y  , 

This  function  was  selected  because  there  is  no  truncation  error,  and  all  of 
the  measured  error  is  due  ;o  inaccuracies  in  the  solution  of  the  difference 
equations.  All  of  the  computations  were  performed  on  a  CDC  66CX!)  computer. 

The  iterative  methods  used  are: 

SOR  :  point  successive  overrelaxation  [25,  p.  58], 

SLOR;  successive  line  overrelaxation  [25,  p.  8o], 

ADI  ;  Peacenjan-Rachford  alternating  direction  implicit  iteration 
[2k ^  Chap,  6], 

The  iteration  parameters  used  are  those  for  the  imbedding  rectangle  R^  ,  and 
for  ADI  the  parameters  for  cycles  of  length  four  are  CEdculated  by  the 
Wachspress  algorithm  [2k ^  Chap.  6],  The  initial  guess  is  identically  zero, 
and  the  iterations  are  terminated  when  the  maximum  difference  between  iterates 
is  less  than  lO"^  , 

The  direct  method  used  is  variant  one  of  the  Buneman  algorithm  [  8  ^  Sec,  ll]. 
Preprocessing  times  are  given  in  Table  2,  Computational  results  for  a  similar 
problem  are  given  in  [1^], 

The  problem  described  in  Section  for  the  region  in  Figure  2  has  been 
treated  by  Daly  and  Nichol3  [10].  The  mesh  used  has  25  x  Ho  =  920  points. 

Using  the  direct  method  of  matrix,  decomposition,  a  particular  solution  requires 
about  50  -  50'5^  of  the  time  required  for  a  point  Gauss-Seidel  iterative  procedure. 


20 


The  prohlen  dlscuesed  in  Section  4  for  the  region  In  Figure  3  hne  been 
treated  by  Horae  and  Budalnakl  [21],  The  »e.h  need  has  52  x  98  .  5096 
polnta,  and  the  preproceaelng  tine  la  approUnately  150  aeconda.  The  region 
and  differential  operator  are  very  aeldon  changed,  ao  the  factored  capacitance 
matrix  la  atored  on  nagnetic  tape.  Thua  there  la  eaaentially  no  preproceaalng 
time  for  the  execution  of  the  progran.  To  aolve  for  a  particular  aolutlon 
requirea  about  2  aeconda,  which  la  approxlBately  Uof,  of  the  time  reoulred  for 
a  successive  line  overrelaxation  iterative  procedure. 


Table  1,  CoBnratational  results  for  solving  the  discrete  Poisson  equation 


Region 


Method 

Maximum 

Error 

Computation 
Time  (Sec.) 

Seeded 

Computation 

Time 

SOR 

5.02  (-6) 

3.586 

21.866 

SLQR 

7.65  (-6) 

2.654 

16.183 

ADI 

2.36  (-6) 

1.128 

6.878 

Direct 

4.44  (-13) 

0.164 

l.OOC 

SOR 

8.12  (-6) 

29.388 

43.994 

SLOR 

7.95  (-6) 

21.424 

52.072 

ADI 

5.41  (-6) 

5.642 

8.446 

Direct 

1.90  (-12) 

0.668 

1.000 

SOR 

2.35  (-6) 

5.570 

21.250 

SLOR 

6.48  (-6) 

2.558 

15.226 

ADI 

2.11  (-6) 

0.870 

5.179 

Direct 

5.77  (-15) 

0.168 

1.000 

SOR 

2.02  (-6) 

29.624 

43.565 

SLOR 

9.96  (-6) 

20.510 

30.162 

ADI 

5.57  (-6) 

5.552 

7.841 

Direct 

1.54  (-12) 

0.680 

1.000 

T*le_2.  Preprocessing  time  for  the  direct 


nethod  results  in  Tehle  1. 


Region  In  the  x-y  plane 


Regl 
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